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I.  INTRODUCTION 


Tunneling  ionization  is  usually  described  using  the  time  dependent  Schrodinger  equation 
(TDSE).  The  TDSE  is  a  strictly  non-relativistic  equation.  In  contrast,  the  dynamics  of  a 
free  particle  become  relativistic  when  the  normalized  vector  potential,  a  =  eA/mc2,  satisfies 
a  >  1.  For  typical  laser  frequencies,  this  corresponds  to  an  irradiance  I  ~  1018  W/cm2. 
Such  irradiances  are  achieved  regularly  in  many  laboratories.  The  corresponding  fields 
are  sufficient  to  strip  nitrogen  down  to  the  K-shell  via  barrier  suppression.  Irradiances 
of  I  >  1020  W/cm2  are  not  achieved  as  frequently,  but  are  easily  within  the  reach  of  several 
ultra-high  power  laser  systems  around  the  world.  It  is  expected  that  these  irradiances  are 
sufficient  to  fully  strip  neon. 

Although  the  dynamics  of  free  electrons  are  relativistic  for  a  >  1,  the  spectra  of  atoms 
that  can  be  brought  to  a  high  charge  state  by  such  a  held  are  non-relativistic.  Put  another 
way,  all  the  bound  electrons  in  an  atom  that  requires  a  >  1  to  be  fully  stripped,  are  well 
described  by  Schrodinger  theory.  To  see  this,  note  that  for  a  hydrogen-like  ion,  the  barrier 
suppression  model  [T]  gives  the  threshold  for  ionization  as 


a  Z3 
UJLTa  16 


(1) 


where  Z  is  the  atomic  number,  is  the  laser  frequency,  ra  =  h3 /me4  is  the  atomic  unit  of 
time,  and  a  =  e2/hc  is  the  fine  structure  constant.  Taking  a  =  1  and  a  laser  wavelength 
A  =  0.8  pm  (cpLTa  =  0.057)  gives  Z  =  5.  Hence,  boron  is  the  heaviest  element  that  is  fully 
stripped  by  a  laser  with  A  =  0.8  pm  and  a  =  1.  The  condition  for  an  atomic  spectrum  to  be 
non-relativistic  is  Z  -C  of  1  «  137,  as  follows  from  elementary  Dirac  theory.  One  concludes 
that  even  though  a  =  1  leads  to  relativistic  motion  of  free  electrons,  the  bound  electrons 
that  are  freed  are  non-relativistic.  This  observation  can  be  useful  when  formulating  the 
initial  conditions  for  a  relativistic  quantum  optics  problem. 

The  process  of  ionization  involves  the  dynamics  of  both  bound  and  free  electrons.  For 
a  >  1  and  Z  <C  137,  the  bound  electrons  are  quantum  mechanical  and  non-relativistic,  while 
the  free  electrons  are  relativistic  and  classical.  If  much  heavier  elements  are  considered,  even 
the  bound  electrons  might  have  to  be  treated  relativistically.  Hence,  a  complete  description 
of  the  problem  must  be  quantum  mechanical  and  fully  relativistic.  In  this  report  we  develop 
a  numerical  tool  capable  of  treating  such  problems. 
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TABLE  I:  Natural  Units 


Quantity 

Unit 

SI  Value 

Velocity 

c 

2.9979  x  108  m/s 

Charge 

\qe  |/a1/2 

1.8755  x  1CT18  C 

Mass 

me 

9.1094  x  10-31  kg 

Length 

h/mec 

3.8616  x  10“13  m 

Time 

h/m.ei :2 

1.2881  x  10~21  s 

Angular  Momentum 

h 

1.0546  x  10-34  J-s 

The  non-relativistic  theory  of  photo-ionization  was  pioneered  by  L.V.  Keldysh  in  1965 
0 .  Shortly  thereafter,  in  the  series  of  papers  by  Perelemov  et  al.  [3]  [6] ,  it  was  brought  into  a 
form  that  has  remained  useful  right  up  to  the  present  day.  Relativistic  photo-ionization  has 
received  less  attention,  although  some  analyses  have  been  given  Fully  time  dependent 

numerical  solutions  of  relativistic  wave  equations  are  only  recently  appearing  in  the  literature 

uni  na¬ 


il.  RELATIVISTIC  WAVE  EQUATIONS 

The  Dirac  equation  describes  the  motion  of  an  electron  in  an  external  potential,  such 
as  the  superposition  of  an  atomic  binding  potential  and  a  laser  field.  Solution  of  the  Dirac 
equation  requires  evolving  a  4-component  bi-spinor  wavefunction.  However,  as  discussed  in 
Ref.  JT2j,  the  Dirac  equation  can  be  separated  into  4  independent  equations  by  neglecting 
terms  involving  spin.  This  results  in  the  Klcin-Gordon  equation 

\(tdt  +  g<F)2  —  (?V  +  qA)2  —  m2]  T  =  0  (2) 

where  q  is  the  charge  and  m  is  the  mass  of  the  particle.  Here,  and  in  all  that  follows,  natural 
units  are  employed  (see  Table  [I]).  Using  the  Coulomb  gauge,  V  •  A  =  0,  and  assuming  a 
static  scalar  potential, 

(□2  -  m2  +  q2AflAfl)  T  +  2 iqd^A^)  =  0  (3) 

where  D2  =  V2  —  <92,  /r  is  a  relativistic  tensor  index,  and  the  metric  signature  is  (H - ). 

Note  that  in  the  chosen  gauge  V  ■  (AT)  =  A  ■  VT.  The  expression  on  the  left  is  useful  for 
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FIG.  1:  Cartoon  of  an  atomic  wavefunction  in  various  scenarios,  (a)  bound  state,  (b)  non- 
relativistic  ionizing  state,  (c)  relativistic  ionizing  state,  (d)  cylindrical  bound  state,  (e)  cylindrical 
non-relativistic  ionizing  state,  (f)  cylindrical  relativistic  ionizing  state.  Coordinate  axes  are  la¬ 
beled  by  basis  vectors  taken  from  any  of  the  Cartesian,  cylindrical,  or  spherical  coordinates,  with 
ignorable  coordinates  in  red. 

finite  volume  differencing,  while  the  one  on  the  right  is  useful  for  analysis. 

Depending  on  the  problem,  various  coordinate  systems  are  used  to  solve  Eq.  (§•  In 
this  report  Cartesian  coordinates  ( x,y,z ),  cylindrical  coordinates  (p,ip,z),  and  spherical 
coordinates  (r,  6,  ip)  are  used.  Cartoons  of  various  wavefunctions  are  shown  in  Fig.  [I]  Panels 
(a),  (b),  and  (c)  are  for  a  central  binding  potential  $(r),  while  (d),  (e),  and  (f)  are  for  a 
cylindrical  binding  potential  $(p).  Naturally,  the  former  case  is  a  more  realistic  model  of  an 
atom  or  ion.  Panels  (a)  and  (d)  show  a  bound  state  wavefunction  for  which  A  =  0.  In  both 
geometries,  two  coordinates  are  ignorable.  Panels  (b)  and  (e)  show  an  ionizing  wavefunction 
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is  the  non-relativistic  dipole  approximation,  where  A  =  A(i)  is  a  function  of  time  only.  The 
effect  of  the  vector  potential  is  to  stretch  the  wavefunction  in  the  polarization  direction.  For 
a  central  potential,  (p  is  ignorable,  while  for  a  cylindrical  potential  z  is  ignorable.  Panels  (c) 
and  (f)  show  an  ionizing  wavefunction  in  the  relativistic  case,  where  the  spatial  dependence 
of  A  leads  to  a  ponderomotive  force  that  stretches  the  wavefunction  in  the  direction  of  the 
photon  momentum,  k.  In  this  case,  there  is  no  coordinate  that  is  ignorable  for  a  central 
potential,  but  there  is  one  that  is  ignorable  for  a  cylindrical  potential.  This  is  the  reason  for 
considering  a  cylindrical  potential. 


III.  STATIONARY  STATES 


A.  Analytical  Solutions  for  Coulomb  Potentials 


The  time  independent  form  of  Eq.  (J3])  is  obtained  by  making  the  substitution  T(r,f) 
0(r)  exp(— tut).  Using  dt§  =  0,  this  gives 


(V2  —  a;2  +  2g<f>u;  +  q2AflAfi  —  m 2)  0  —  2 iqS7  ■  (A 0)  =  0  (4) 


Specializing  to  the  case  of  a  uniform  magnetic  held  A  =  ^pB^e^  gives 

^V2  —  a;2  +  2q$u  +  g2$2  +  £zqB0  —  i p2q2Bl  —  m2^  R(p,  z)  —  0  (5) 

Here,  the  cylindrical  coordinates  are  denoted  (p,<p,z),  and  iz  is  the  magnetic  quantum 
number,  defined  such  that 

0(r)  =  R(p,  z )  exp(y4 tp)  (6) 


As  discussed  above,  the  use  of  cylindrical  coordinates  is  motivated  by  the  fact  that  for  a 
cylindrical  atom,  one  coordinate  is  ignorable  even  in  the  case  of  a  fully  relativistic  problem. 
The  effect  of  this  fictitious  geometry  on  the  energy  levels  is  determined  below. 

In  the  case  of  a  Coulomb  potential  with  B0  =  0,  the  bound  states  can  be  determined 
exactly.  Approximate  solutions  can  be  found  for  a  weak  magnetic  held.  Consider  hrst 
a  spherical  potential,  $  =  Q/r,  where  Q  =  Za 1//2  is  the  charge  of  the  nucleus  (we  are 
considering  hydrogen-like  ions).  Then  the  bound  state  energies  are 


uj  —  m 


-1/2 


1  + 


QV _ ^ 

V^  +  i)2- w)2 


(7) 
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FIG.  2:  Energy  levels  of  spherical  ions  according  to  Schrodinger,  Klein-Gordon,  and  Dirac  theories, 
for  (left)  hydrogen-like  neon  and  (right)  hydrogen-like  erbium.  The  magnetic  field  is  zero. 


where  c oc  =  qB^/m  is  a  signed  cyclotron  frequency,  nr  is  the  radial  quantum  number  and  £ 
is  the  orbital  quantum  number.  The  principle  quantum  number  is  n  =  nr  +  £  +  1,  which 
completely  determines  the  energy  in  the  non-relativistic  limit,  if  B0  =  0.  The  radial  eigen¬ 
functions,  to  within  a  normalization  factor,  are 


R(r )  =  F{-nr,  2(a  +  1),  2 kr)rae~kr 


(8) 


where  F  is  the  confluent  hypergeometric  function,  a  =  — 1/2  ±  yf{£  +  1/2)2  —  Q2q 2,  and 
k  =  \/m2  —  £zmu>c  —  a;2.  The  condition  that  the  field  should  be  weak  reads  y/uJ^  -C  k. 

A  comparison  of  the  energy  levels  predicted  by  Klein-Gordon  theory,  Schrodinger  theory, 
and  Dirac  theory,  is  shown  in  Fig.  [2j  The  levels  shown  are  taken  from  the  set  {nr  G 
0, 1,  2,  3, 4}  G  0, 1, 2,  3, 4}.  In  the  Schrodinger  theory  the  energy  depends  only  on  the 
combination  nr  +  t.  Note  that  for  neon,  the  difference  is  subtle,  whereas  for  erbium  it  is 
pronounced.  For  elements  heavier  than  erbium,  the  square  root  in  the  expression  for  u 
becomes  imaginary,  so  that  no  stationary  state  exists  for  £  =  0  and  Z  >  68. 

For  a  cylindrical  potential,  $  =  Q/p  [T6],  the  bound  state  energies  are 

r  1  -V2 
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FIG.  3:  Energy  levels  according  to  spherical  Schrodinger  theory  (blue  dashed),  spherical  Klein- 
Gordon  theory  (red),  and  cylindrical  Klein-Gordon  theory  (green)  for  hydrogen-like  neon  with 
B0  =  0. 

and  the  eigenfunctions  are 

R(p)  =  F(-nr ,  2(6  +  1),  2 kp)pbe~kp  (10) 

where  6  =  —  Q2q2.  Evidently  for  cylindrical  hydrogenic  ions,  there  are  no  solutions 

with  tz  =  0  for  any  Z .  A  comparison  of  the  energy  levels  for  cylindrical  and  spherical  neon 
is  shown  in  Fig.  [3] 


B.  Numerical  Solutions  for  Soft  Core  Potentials 


In  numerical  problems,  the  Coulomb  potential  can  only  be  approximated  due  the  singular¬ 
ity  at  the  origin.  Typically  one  replaces  the  Coulomb  potential  with  a  “soft-core  potential,” 
which  for  cylindrical  atoms  has  the  form 


<f> 


Q 

\J  5p2  +  p 2 


(11) 


Obviously,  as  5p  — >  0  the  soft-core  potential  becomes  a  Coulomb  potential.  In  order  to  solve 
the  nonlinear  eigenvalue  problem  (|5|  with  a  soft-core  potential,  a  numerical  scheme  has  to 
be  employed.  Due  to  the  quadratic  form  of  the  nonlinearity,  the  system  can  be  reduced  to 
a  linear  problem  of  twice  the  size  H3-  In  particular,  let  R(p)  be  discretized  on  a  sequence 
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of  mesh  points  {p;  =  (2 i  —  l)4r| *  G  N}.  Evaluating  the  Laplacian  operator  in  (JHj)  by  finite 
volumes  gives  the  block  tridiagonal  matrix  equation 


(12) 


where  R  =  (Rq,  Ri,  R2,  ■  ■ .),  Ri  =  R(pi),  R'  =  ujR,  T  is  a  tridiagonal  matrix,  and  D  is  a 
diagonal  matrix.  The  non-zero  elements  of  T  and  D  are 


Ti,i- 1  — 


Pi  -  Ap/2 

pAp2 


(13a) 


T  = 

- L  2,2 


Pf  Aps 


+  4gA 


(13b) 


_  pi  +  Ap/2 

pAp2 


(13c) 


A,<  =  (13d) 

where  <4  =  ‘h(pi).  Standard  sparse  matrix  packages  can  solve  this  system  for  various  subsets 
of  the  eigenvalue-eigenvector  pairs.  A  similar  system  can  be  constructed  for  spherical  atoms. 

One  important  difference  between  the  bound  states  of  a  soft-core  potential  compared 
with  those  of  a  Coulomb  potential  is  that  a  solution  exists  for  tz  —  0  in  the  former  case, 
but  not  the  latter.  This  is  illustrated  in  Fig.  [4j  which  plots  the  nr  =  0  energy  for  tz  —  0 
and  4  =  1  as  a  function  of  the  soft-core  radius,  bp.  The  4  =  1  energy  is  insensitive  to  bp, 
while  the  4  =  0  energy  exhibits  singular  behavior  as  bp  — »  0.  Of  course,  the  energy  of  the 
discretized  system  does  not  actually  diverge,  even  when  bp  =  0,  due  to  the  fact  that  the 
grid  is  constructed  so  that  the  point  p  =  0  is  never  sampled. 

A  discussion  of  divergent  energies  in  a  Coulomb  potential  can  be  found  in  Ref.  [15]. 


C.  Scalar  Zeeman  Effect 


Numerical  solution  of  the  nonlinear  eigenvalue  problem  (12)  allows  one  to  characterize  the 


splitting  of  the  energy  levels  of  high-Z  ions  in  extreme  magnetic  fields  for  scalar  wavefunc- 
tions.  Consider  hydrogen-like  cylindrical  ununoctium  (Z=118,  the  heaviest  known  element). 
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FIG.  4:  Two  lowest  cylindrical  bound  state  energies  with  Z  =  54  as  a  function  of  soft  core 
radius,  dp.  Dashed  curve  is  cuqo  =  u(nr  =  0,£z  =  0)  and  solid  curve  is  ojq \ .  The  quantity 


koi  =  \Jm?c2 /h2  —  uj^/c2  characterizes  the  size  of  the  ion  (see  Eq.  10). 


FIG.  5:  Weak  field  (dashed)  and  numerical  (solid)  solution  for  energy  levels  in  hydrogen-like 
ununoctium  ( Z  =  118)  with  £z  =  1  (upper  branches)  and  £z  =  —1  (lower  branches). 

For  the  numerical  calculation,  the  number  of  mesh  points  used  is  213,  while  the  mesh  spacing 
and  soft-core  radius  are  ko5p  =  koAp  =  0.00125,  where  ko  =  k{ujc  =  0).  A  comparison  of 
the  numerical  solution  and  the  weak-field  solution  is  shown  in  Fig.  [5]  The  characteristic  size 
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pmc/h 


FIG.  6:  Radial  eigenfunctions  in  a  strong  field  ( Bq  =  10)  for  nr  =  0  (red),  nr  =  2  (blue),  and 
nr  =  4  (green),  with  Z  =  118  and  iz  =  1. 


of  the  ununoctium  ion  is  1/fc  ~  1.5,  so  that  the  weak  held  limit  corresponds  to  ujc  <C  0.44. 
Inspection  of  Fig.  [5]  shows  that  if  uc  is  an  order  of  magnitude  below  the  weak  held  cutoff, 
the  two  solutions  agree.  At  larger  helds,  the  weak  held  solution  underestimates  the  shift  for 
the  upper  branch,  and  overestimates  it  for  the  lower  branch.  A  few  radial  eigenfunctions  in 
a  strong  held  are  shown  in  Fig.  [6j 

In  ordinary  units,  the  weak  held  cutoff  of  ujc  =  0.44  corresponds  to  Bq  «  20  TG,  with 
the  peculiar  result  that  Bq  =  1  TG  is  a  weak  held.  At  present,  such  helds  are  observed  only 
in  connection  with  astrophysical  phenomena  [Qj.  The  highest  laboratory  helds  to  date  are 
laser  generated,  and  are  still  limited  to  hundreds  of  MG  021.  Our  primary  interest  here  is  in 
code  validation,  i.e.,  establishing  results  that  can  be  compared  with  the  fully  time-dependent 
calculations  discussed  below. 

In  the  strong  held  limit,  one  can  make  progress  analytically  by  means  of  the  quasi-classical 
approximation.  Defining  u(p)  =  R(p)p1^2,  the  eigenvalue  problem  takes  the  form  of  a  one 
dimensional  Schrodinger  equation 


u"  +  2  [e(u;)  —  U (a p)\u  =  0 

where  the  effective  energy  eigenvalue  is  e(pj)  =  tn2/2  —  izmu>c/2 
potential  is 


U  (u,  p) 


uQq  b 2 


1/4 


2  „2 


+  g  W 


(14) 


m2 / 2  and  the  effective 


P 


2  p1 


(15) 
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FIG.  7:  Comparison  of  Bohr-Sommerfeld  energy  levels  (circles)  with  solution  of  discretized  eigen- 
system  (squares),  for  Bq  =  0.5,  £z  =  1,  and  Z  =  118. 


For  states  with  nr  1  or  tz  1  (equivalently,  uj  — >  m),  Bohr-Sommerfeld  quantization 
may  be  employed.  The  quantization  rule  is 


'pi 


P2  _  /  i ' 

<W2(eM  -UifZp) 1  =  I  nB  +  j  I  * 


(16) 


where  ns  +  1  G  N,  and  p±  and  p2  are  roots  of  the  integrand.  We  associate  the  Bohr- 
Sommerfeld  quantum  number,  ub,  with  the  radial  quantum  number,  nr.  Fig.  [7]  compares 
the  energy  levels  for  B0  =  0.5,  £z  =  1,  and  Z  =  118,  as  computed  using  Bohr-Sommerfeld 


quantization  and  direct  numerical  solution  of  the  eigensystem  (12).  As  expected,  the  agree¬ 
ment  improves  for  higher  quantum  numbers.  The  error  asymptotes  to  a  finite  value  due  to 
discretization  errors  in  the  solution  of  (12).  Fig.  [8]  shows  a  similar  comparison  for  several 
values  of  the  magnetic  field,  holding  ub  fixed  at  ns  =  10.  For  the  highest  value  of  the 
magnetic  field,  the  solutions  for  several  values  of  iz  are  displayed. 


IV.  TIME  DEPENDENT  SIMULATIONS 

A.  Numerical  Algorithm 

The  Klein-Gordon  equation  is  a  second  order  hyperbolic  equation.  It  can  be  explicitly 
differenced,  provided  a  Courant-type  condition  is  satisfied.  Centered  time  differencing  of 
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FIG.  8:  Comparison  of  Bohr-Sommerfeld  energy  levels  (circles)  with  solution  of  discretized  eigen- 
system  (plus  signs)  as  a  function  of  magnetic  field,  for  ng  =  10,  £z  =  1,  and  Z  =  118.  For  the 
highest  value  of  uc,  solutions  are  plotted  for  —  3  <  iz  <  3 


Eq.  (j3j),  with  <9f<3>  =  0,  gives 


Vj/n+i 

At2 


(V2  +  g2<f>2  —  q2An 2  —  m2)\&n  —  2 zV  ■  (gAra\I/n)  2Tn  —  (1  —  iq&At)'# 


n—  1 


1  +  iq&At 


+ 


At2(l  +  iq&At) 


(U) 


where  n  indexes  the  time  levels.  The  spatial  derivatives  can  be  put  in  finite  difference  form 
by  means  of  finite  volumes.  On  a  two-dimensional  Cartesian  grid,  this  results  in 


hj 


T”  ,  ■ 

*- ij 


-  2Tn  ■ 

h3 


+ 


i+lp 


Ax 2 


+ 


i  — 2Tn  • 


+ 


*ji+l 


Ay2 


(18a) 


v  ■  *  a;1  t;'  )  = 


a,  ,.,>1/''  A,-.,-.,*;:,..,  -  A;„  ,vi/;v  , 


+ 


(18b) 


2Ax  2A  y 

where  i  and  j  index  the  mesh  points.  The  extension  to  three  dimensions  is  straightforward. 
The  numerical  properties  of  this  scheme  can  be  evaluated  by  inserting  the  form  U /”  •  = 

exp[z(ipxAx  +  jpyAy  +  kpzAz  —  nujAt)]  into  the  difference  equations,  with  a  constant, 
uniform  A  and  $.  The  resulting  numerical  dispersion  relation  is 


sin2  =  F( p)  +  G(w)  +  At2  (l n 2  +  q2A2  -  qV) 


(19) 
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where 

f  (p)  =  £ 

i 

G(u)  =  -Atq$  sin  (a;  At)  (20b) 

In  the  equation  for  F{ p),  i  indexes  Cartesian  coordinates,  with  Ai  =  Ax,  etc..  It  is  easily 
verified  that  in  the  limit  where  the  arguments  of  the  trigonometric  functions  are  small,  the 
numerical  dispersion  relation  reduces  to  the  true  dispersion  relation, 

(a;  —  g<P)2  =  (p  —  qA )2  +  m2  (21) 

For  a  free  particle,  A  =  $  =  0,  and  the  energy  will  always  be  real  provided 

1/1  1  1 
At 2  \Ax2  Ay 2  +  A z2 

In  one  dimension,  this  becomes  At  <  Az(l  +  m2 Az2 / A)~1/'2 .  Note  that  if  the  cell  size  is 
much  smaller  than  the  Compton  wavelength  of  the  particle  in  question  (Az  <C  1/m),  then 
the  familiar  Courant  condition  for  the  wave  equation  is  recovered.  A  more  general  stability 
criterion  is  0  <  Ai  <  1,  where 

M  =  F( p)  +  G(u>)  +  ^At2(m2  +  q2A2  -  g2<P2)  (23) 

Here,  unlike  the  case  of  a  free  particle,  instabilities  can  occur  either  in  the  large  or  small 
momentum  limits.  This  is  because  even  when  p  — >  0,  a  large  <3>  can  still  drive  M.  negative. 
The  one-dimensional  numerical  dispersion  relation  for  g<3>  =  —10  and  qA  =  0  is  shown  in 
Fig.  0  Discretization  parameters  leading  to  instability  give  Fig.  |9](a),  while  parameters 
leading  to  stability  give  Fig.  [9](b). 

B.  Validation  Against  Zeeman  Effect 

Two  codes  have  been  written  that  implement  the  above  algorithm  on  General  Purpose 
Graphical  Processing  Units  (GPGPU).  One  is  programmed  using  Python  and  OpenCL,  and 
the  other  is  a  turboWAVE  module  programmed  using  C++  and  OpenCL.  The  turboWAVE 
version  supports  combined  Message  Passing  Interface  (MPI)  and  GPGPU  programming.  It 
is  fully  three  dimensional,  although  the  examples  treated  here  are  two-dimensional. 


At2 
A J 


At2 

sin2 (p.iAi/2)  -  —qAiSmfaAi) 


(20a) 
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(a)  (b) 


FIG.  9:  Numerical  vs.  actual  dispersion  relation  with  (/<!>  =  —10  and  qA  =  0.  Blue  curves  are 
the  true  energy,  red  squares  are  the  real  part  of  the  numerical  energy,  and  green  curves  are  the 
imaginary  part  of  numerical  energy.  Panel  (a)  is  for  Az  =  0.2  and  At  =  0.1,  while  panel  (b)  is  for 
Az  =  0.02  and  At  =  0.01. 


In  order  to  validate  the  two  codes,  they  are  benchmarked  against  time  independent 
calculations  of  relativistic  Zeeman  splitting.  In  order  to  extract  the  spectrum  of  stationary 
states  from  a  time  dependent  code,  it  is  desirable  to  use  an  initial  condition  that  is  a 
superposition  of  all  possible  states.  At  the  same  time,  the  results  should  not  be  prejudiced 
by  the  initial  condition.  An  entirely  random  initial  wavefunction  serves  both  purposes.  After 
allowing  this  wavefunction  to  evolve  for  a  suitably  long  time,  the  spectrum  can  be  estimated 
by  Fourier  transformation  in  time  at  selected  points  in  space.  The  spectrum  at  each  spatial 
point  is  different  in  terms  of  the  relative  magnitude  of  the  spectral  lines,  but  the  energy  of 
each  line  is  independent  of  the  spatial  point  chosen. 

Two  simulations  are  performed,  one  with  Bq  =  0  and  one  with  Bq  =  0.5.  The  binding 
potential  is  a  cylindrical  soft-core  potential,  with  Z  =  118  and  5p  =  0.2.  The  grid  contains 
8000  x  8000  cells,  with  cell  size  0.02  x  0.02.  The  time  step  is  0.01  and  9  x  105  steps  are 
taken.  The  results  for  Bq  =  0  and  Bq  =  0.5  are  shown  in  Fig.  [lOj  The  solid  curves  in 
the  form  of  discrete  peaks  are  the  Fourier  transformed  wavefunction  at  two  different  spatial 
points.  As  expected,  the  height  of  the  peaks  depends  on  the  spatial  point,  but  the  energy 


does  not.  The  energies  are  also  determined  by  solving  the  time  independent  system  (12)  for 


various  values  of  nr  and  tz.  These  solutions  are  shown  as  vertical  dashed  lines.  Except  for 
the  highest  lying  bound  states  in  the  Bq  =  0  case,  the  two  calculations  agree.  Actually,  a 
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FIG.  10:  Comparison  of  energy  levels  in  a  cylindrical  soft  core  potential  with  Z  =  118,  with  (a) 
Bq  =  0,  and  (b)  Bq  =  0.5,  as  determined  by  time  dependent  (solid)  and  time  independent  (dashed) 
calculations.  The  brown  lines  are  the  Fourier  spectrum  evaluated  at  ro  =  2ex  and  the  blue  lines 
are  the  Fourier  spectrum  evaluated  at  ro  =  5ex.  The  blue  spectrum  is  multiplied  by  four. 

real  ion  has  an  unbounded  density  of  states  as  u  — >  1  when  B0  =  0,  whereas  any  numerical 
ion  has  only  a  finite  number  of  states.  It  should  also  be  noted  that  in  a  magnetic  held,  the 
bound  state  spectrum  is  not  bounded  from  above,  i.e. ,  bound  states  exist  for  u  >  1. 
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TABLE  II:  Parameters  for  simulation  of  relativistic  ionization 


Parameter 

Symbol 

Value 

Steps 

Nt 

200000 

Time  Step 

At 

0.2  h/mc2 

Cells 

Nx  X  Ny 

4000  x  12000 

Space  Step 

Ax,  Ay 

h/mc 

Residual  Charge 

Q 

18a1/2 

Soft  Core  Radius 

5p 

h/mc 

Radial  Quantum  Number 

nr 

0 

Magnetic  Quantum  Number 

4 

1 

Laser  wavelength 

A 

24  nm 

Vector  Potential 

do 

3.3 

C.  Relativistic  Ionization  Example 

Relativistic  photo-ionization  experiments  will  likely  involve  laser  radiation  at  a  wave¬ 
length  of  A  =  0.8  or  A  =  1.05  /jm.  In  order  to  model  one  cycle  of  the  radiation,  the  number 
of  time  steps  has  to  satisfy  Nt  A/Ac  ~  105,  where  Ac  is  the  Compton  wavelength.  The 
number  of  grid  cells  in  one  dimension  has  to  satisfy  the  same  inequality  due  to  the  fact  that 
the  spatial  scale  of  the  orbit  of  a  classical  electron  in  a  plane  wave  is  A  in  the  relativistic 
limit.  Based  on  the  performance  characteristics  given  below,  a  single  GPGPU  would  take 
about  106  hours  to  complete  such  a  calculation,  even  in  two  dimensions  (assuming  it  had 
large  enough  memory).  As  a  result,  thousands  of  GPGPU  compute  nodes  would  be  required 
to  simulate  a  relativistic  photo-ionization  experiment  at  full  scale. 

In  order  to  run  an  example  on  the  currently  available  GPGPU  cluster  “Dirac”  at  the 
National  Energy  Research  Supercomputing  Center  (NERSC),  the  parameters  displayed  in 
Table  [U]  are  used.  The  simulation  runs  in  about  1.2  hours.  The  primary  scale  reduction  is 
the  fictitious  laser  wavelength  of  24  nm.  The  increased  laser  frequency  means  that  fewer 
steps  and  grid  cells  have  to  be  used  to  model  a  full  cycle  of  the  electron  motion.  The  laser 
field  has  the  form 


A (z,  t )  =  a0  [cos(u;i>  —  —  1]  0  (t  —  z)  ey 


(24) 
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where  0  is  the  Heaviside  step  function.  The  electrostatic  field  is  a  soft  core  potential.  The 
value  of  a0  =  3.3  is  determined  by  calculating  the  barrier  suppression  threshold  [T]  for  an  ion 
with  residual  charge  18a1/2  (cylindrical  H-like  argon)  and  ground  state  energy  0.9961.  The 
ground  state  energy  is  determined  from  Eq.  with  nr  =  0  and  tz  —  1,  and  corresponds 
to  an  ionization  potential  of  about  2  keV.  An  interesting  question  is  how  well  the  barrier 
suppression  model,  which  is  based  on  electrostatic  arguments,  works  in  the  relativistic  limit. 

The  simulated  ionization  rate  is  defined  in  terms  of  the  charge  current  flowing  out  of 
a  volume  containing  the  bound  state.  The  Klein-Gordon  equation  satisfies  the  continuity 
equation  dfg  +  V  ■  j  =  0,  where 

g  =  —  |T|2  (25) 

Tfl 

is  the  analog  of  the  classical  charge  density  and 

2 

j  =  — —  (T*VT  -  TVT*)  -  —  ITI2  (26) 

2  m  y  m 

is  the  analog  of  the  classical  current  density.  The  expectation  value  of  the  charge  contained 

in  a  ball  B  is 

(q)b=  [  d3re  (27) 

Jb 

Take  the  radius  of  B  to  be  large  enough  so  that  it  contains  the  charge  associated  with  the 
bound  state  almost  entirely.  Then  the  ionization  probability  is  defined  as  the  expectation 
value  of  the  charge  outside  B ,  divided  by  the  total  charge  of  the  particle: 


P  = 


1  -  (l)i 


(28) 


The  ionization  rate  is 


dt  q  dt 


d3rg 


(29) 


which  by  the  divergence  theorem,  is  the  same  as  the  current  flowing  out  of  the  volume, 
divided  by  the  charge.  It  should  be  noted  that,  due  to  the  fact  that  g  can  take  either  sign, 
this  definition  can  lead  to  negative  ionization  probabilities.  If  the  system  being  considered 
is  an  ordinary  ion,  one  may  interpret  P  <  0  as  an  indication  that  finding  a  positron  outside 
of  B  is  more  likely  than  finding  an  electron.  An  intriguing  possibility  is  that  the  laser  held 
produces  a  free  positron  and  an  additional  bound  electron. 

The  ionization  rate  as  defined  above  is  displayed  in  Fig.  [TT|)a).  The  applied  electric  held 
is  displayed  in  Fig.  [TT|(b) .  Some  of  the  characteristics  of  the  ionization  rate  are  familiar  from 
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(a) 


Time  (as) 


FIG.  11:  Simulation  of  relativistic  photo-ionization. 


(a)  Rate  defined  as  in  (29)  where  B  is  a 


cylinder  of  radius  2.3  A,  (b)  Applied  electric  field. 


(a) 


FIG.  12:  Charge  density  associated  with  ionizing  relativistic  wavefunction  evaluated  at  (a)  t  =  26 
as  (b)  t  =  39  as  (c)  t  =  52  as. 


the  non- relativistic  tunneling  theory.  First,  the  total  ionization  probability  (area  under  the 
rate  curve)  is  in  the  range  of  a  few  percent,  so  that  the  threshold  field  from  [1]  does  indeed 
serve  as  a  suitable  threshold  for  ionization.  Moreover,  the  peak  of  the  ionization  rate  occurs 
near  the  peak  of  the  electric  field,  as  ordinary  tunneling  theory  predicts.  An  unexpected 
feature  is  that  the  ionization  rate  has  multiple  peaks  within  one  half-cycle  of  the  applied 
field. 
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In  order  to  visualize  the  relativistic  ionizing  wavefunction,  the  charge  density,  g,  is  plotted 


as  a  false-color  image  in  Fig.  12,  The  three  panels  show  g  evaluated  at  three  different  times. 
The  three  peaks  visible  in  the  ionization  rate  correspond  to  the  white  and  blue  islands  visible 


against  the  dark  background  in  Figs.  12  b)  and  (c).  The  remnant  of  the  bound  state  is  the 
dark  red  feature  at  (x,y)  =  (0,0).  One  can  clearly  see  a  feature  of  the  Volkov  solution 
for  an  electron  in  a  plane  wave.  Namely,  the  wavefunction  is  bent  in  the  direction  of  the 
electromagnetic  wave  propagation.  One  would  like  to  characterize  how  this  “ponderomotive” 
effect  alters  the  recollision  phenomena  that  are  well  known  in  the  non-relativistic  case. 


V.  COMPUTATIONAL  PERFORMANCE 


As  discussed  above,  the  fundamental  scale  ratio  in  a  relativistic  quantum  optics  problem 
is  A/Ac,  where  A  is  the  laser  wavelength  and  Ac  is  the  Compton  wavelength.  Due  to  the 
large  order  of  magnitude  of  this  ratio,  computational  performance  is  crucial.  To  this  end, 
we  developed  a  relativistic  quantum  optics  module  within  the  turboWAVE  framework  m- 
Among  other  things,  turboWAVE  provides  a  programming  environment  for  solving  partial 
differential  equations  on  various  structured  grids.  The  framework  provides  abstractions  that 
are  useful  for  parallelization.  Parallelization  across  distributed  processors  is  accomplished 
via  the  Message  Passing  Interface  (MPI).  Parallelization  across  shared  memory  processors 
is  accomplished  via  the  OpenCL  language.  In  the  latter  case,  particular  attention  is  given 
to  general  purpose  graphical  processing  units  (GPGPU). 

Parallelization  of  the  Klein-Gordon  equation  is  accomplished  via  domain  decomposition, 
where  each  GPGPU  advances  the  solution  in  a  given  domain,  and  MPI  is  used  for  commu¬ 
nication  between  domains.  The  structure  of  a  two-dimensional  domain  is  shown  in  Fig.  13 
(the  extension  to  three  dimensions  is  straightforward,  and  has  been  implemented  in  tur¬ 
boWAVE).  Cells  are  labeled  either  by  an  index  space  pair,  which  correspond  to  a  spatial 
location,  or  a  single  index  that  corresponds  to  a  location  in  computer  memory.  The  region 
of  memory  where  the  domain  resides  is  called  the  compute  buffer.  In  the  figure,  the  white 
layer  of  cells  are  ghost  cells,  the  blue  cells  are  edge  cells,  and  the  orange  cells  are  interior 
cells.  The  ghost  cells  and  edge  cells  together  are  called  boundary  cells.  The  interior  and 
edge  cells  can  be  updated  independently  on  each  domain  provided  the  ghost  cells  contain 
the  values  of  the  edge  cells  in  adjacent  domains.  The  problem  of  parallelizing  the  solution, 


19 


30 

31 

32 

33 

34 

35 

24 

25 

26 

27 

28 

29 

18 

19 

20 

21 

22 

23 

12 

13 

14 

15 

16 

17 

6 

7 

8 

9 

10 

11 

0 

1 

2 

3 

4 

5 

0  1  2  3  4  5 


Index  Space  Axis  0 


FIG.  13:  Structure  of  a  typical  computational  domain.  The  index  space  axes  correspond  directly 
to  physical  axes  in  space,  while  the  indices  written  in  each  cell  correspond  to  a  location  in  memory. 

then,  reduces  to  that  of  updating  the  ghost  cells  in  the  memory  of  a  given  GPGPU  using 
information  from  other  GPGPUs.  An  efficient  realization  of  this  operation  is  outlined  as 
follows: 

1.  Device:  Copy  boundary  cells  from  compute  buffer  into  a  contiguous  transfer  buffer 

2.  Device:  Copy  transfer  buffer  to  host  memory 

3.  Host:  Update  ghost  cells  using  MPI  as  usual 

4.  Device:  Read  updated  transfer  buffer  from  host  memory 

5.  Device:  Copy  cells  from  transfer  buffer  to  compute  buffer 

Here,  the  OpenCL  terminology  “device”  and  “host,”  typically  refers  to  a  GPGPU  and  CPU, 
respectively.  Operations  with  the  transfer  buffer  are  facilitated  by  storing  an  index  map  on 
the  GPGPU.  The  index  map  contains  pairs  of  indices,  where  one  index  in  the  pair  points  to 
a  location  in  the  compute  buffer,  and  the  other  points  to  the  corresponding  location  in  the 
transfer  buffer. 

Once  the  ghost  cells  have  been  updated,  the  GPGPU  can  advance  the  relativistic  wave- 
function  using  the  Klein-Gordon  equation.  This  requires  that  the  compute  buffer  contain 
Tn_1,  Tra,  and  A™  [19].  In  cases  where  the  potential  is  time  varying,  an  efficient  way  of 
updating  it  is  required.  In  order  to  avoid  calling  special  functions  on  the  GPGPU,  both  the 
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FIG.  14:  Performance  in  terms  of  millions  of  cells  (Mcells)  per  second  for  various  (a)  CPU  cores 
and  (b)  GPGPUs.  Blue  bars  are  for  single  precision  and  green  bars  are  for  double  precision. 

Python  code  and  turbo  WAVE  assume  A  is  in  the  form  of  a  plane  wave.  The  host  updates 
a  single  longitudinal  strip,  which  is  then  copied  to  a  device  buffer.  The  GPGPU  then  dupli¬ 
cates  the  strip  at  every  transverse  index.  The  Klein-Gordon  update,  as  described  in  section 
TVA  is  then  executed.  Each  cell  in  the  compute  buffer  is  an  OpenCL  work  item.  Since  the 


algorithm  involves  no  spatial  derivatives  of  it  can  be  carried  out  in  place  by  storing 

Tn+1  in  the  memory  occupied  by 

The  performance  of  the  algorithm  on  a  single  compute  node  (no  network  message  passing) 


is  illustrated  in  Fig.  14  for  several  CPU  and  GPGPU  devices.  One  notable  characteristic  is 
that  the  single  and  double  precision  performance  is  nearly  the  same  on  a  CPU,  but  differs 
by  more  than  two-fold  on  a  GPGPU.  In  terms  of  double  precision  performance,  a  GPGPU 
is  typically  100  times  faster  than  a  single  CPU  core.  In  the  most  extreme  case,  the  NVIDIA 
GTX  680  running  in  single  precision  mode  is  over  700  times  faster  than  one  core  of  the 
Opteron  6172  “Magny- Corns.”  The  floating  point  efficiency  (sustained  to  peak  throughput) 
ranges  between  4%  and  10%  for  both  CPU  cores  and  GPGPUs. 


Finally,  Fig.  [15]  shows  the  results  of  a  scaling  study  carried  out  on  the  GPGPU  cluster 
“Dirac”  at  the  National  Energy  Research  Supercomputing  Center  (NERSC).  The  study  was 
carried  out  up  to  the  maximum  number  of  GPGPU  nodes  that  are  made  available  by  the 
queuing  system.  Fig.  [I5|(a)  shows  the  performance  in  the  aggregate  vs.  the  number  of  MPI 
compute  nodes,  where  each  node  is  assigned  a  single  GPGPU.  The  results  are  shown  for 
three  different  size  grids.  The  results  normalized  to  the  number  of  GPGPUs  are  shown  in 
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FIG.  15:  Scaling  study  on  NERSC  GPGPU  cluster  “Dirac”  for  three  different  size  grids  showing  (a) 
millions  of  cells  per  second  in  the  aggregate,  and  (b)  millions  of  cells  per  second  on  each  GPGPU. 
Double  precision  is  used  in  all  cases.  One  GPGPU  is  assigned  to  each  MPI  node. 


Fig.  15  b).  The  performance  per  GPGPU  on  the  smallest  grid  drops  noticeably  in  going  from 
1  to  12  nodes.  For  the  larger  size  grids,  the  performance  is  well  sustained.  An  interesting 
feature  is  that  the  performance  per  GPGPU  actually  increases  in  going  from  1  to  2  nodes 
on  the  medium  grid.  This  may  indicate  there  is  a  performance  penalty  incurred  when  the 
device  memory  is  nearly  exhausted.  Note  that  for  the  large  grid,  the  case  of  a  single  node 
could  not  be  carried  out  due  to  reaching  the  device  memory  limit.  In  the  best  case  each 
GPGPU  achieves  45  GFLOPS  per  second,  about  9%  of  the  theoretical  maximum.  The  best 
floating  point  performance  in  the  aggregate  is  0.5  TFLOPS  per  second. 


VI.  CONCLUSIONS 

The  availability  of  high  performance  computing  resources  makes  it  possible  to  carry  out 
relativistic  quantum  optics  calculations  from  first  principles.  A  starting  point  is  to  solve  the 
Klein-Gordon  equation  for  a  spin  zero  particle.  Solutions  of  the  time  independent  equation 
are  developed  analytically  and  numerically  in  order  to  initialize  time  dependent  calcula¬ 
tions.  Energy  levels  extracted  from  the  time  dependent  code  are  in  agreement  with  those 
obtained  from  time  independent  calculations  and  analytical  solutions.  The  time  dependent 
calculation  provides  solutions  to  relativistic  quantum  optics  problems  such  as  relativistic 
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photo-ionization  of  high-Z  ions  by  extreme  fields.  This  type  of  calculation  is  made  possi¬ 
ble  by  modern  heterogeneous  computing  resources,  and  opens  up  the  little-explored  field  of 
relativistic  quantum  optics. 
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